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Abstract 



We study pattern formation processes in anisotropic system governed by the Kuramoto-Sivashinsky equation with 
multiplicative noise as a generalization of the Bradley-Harper model for ripple formation induced by ion bombard- 
ment. For both linear and nonlinear systems we study noise induced effects at ripple formation and discuss scaling 
behavior of the surface growth and roughness characteristics. It was found that the secondary parameters of the 
ion beam (beam profile and variations of an incidence angle) can crucially change the topology of patterns and the 

ryj corresponding dynamics. 

+J PACS 05.40.-a, 05.70.Ln, 64.60.1, 68.35. Ct, 79.20.Rf 

i 1 Introduction 

A fabrication of nanoscale surface structures have attracted a considerable attention due to their applications 
in electronics pQ. In the last five decades many studies have been devoted to understanding the mechanism of 

l — ' pattern formation and its control during ion-beam sputtering (see, for example Refs . [21 El HI El El EZl [8 , 9 ]). Among 
theoretical investigations there are a lot of experimental data manifesting a large class of patterns appeared as 
result of self-organization process on the surface of a solid. It was shown experimentally that main properties of 

(■■<") pattern formation and structure of patterns depend on the energetic ion-beam parameters such as ion flux, energy 

of deposition, angle of incidence and temperature. Formation of ripples was investigated on different substrates, i.e. 
on metals (Ag and Cu) DUE] on semiconductors (Ge pLj and Si [331 HH [15] ) on Sn [16] , InP [32], on Cd 2 Nb 2 7 
pyrochlore |18] and other. Height modulations on the surface induced by ion-beam sputtering result in formation 
of ripples having the typical size of 0.1 to 1 [im and nanoscale patterns with the linear size of 35 to 250 A [15] . 

It is well known that orientation of ripples depends on the incidence angle. At the incidence angles around 
7r/2 the wave-vector of the modulations is parallel to the component of the ion beam in the surface plane, whereas 
at small incidence angles (close to grazing) the wave- vector is perpendicular to this component. The orientation 
of ripples can be controlled by a penetration depth which is proportional to the deposited energy. Analytical 
investigations provided by Cuerno and Barabasi show a possible control of pattern formation governed by both the 
incidence angle and penetration depth [H [5] . The main theoretical models describing ripple formation are based 
on results of the famous works of Bradley and Harper [3 , Kardar, Parisi, and Zhang [20], Wolf and Villian [2"T] . 
Kuramoto and Sivashinsky j22j . The main mechanisms for pattern formation were set to predict orientation change 
of the ripples, formation of holes and dots. These models were generalized taking into account additive fluctuations 
leading to statistical description of the corresponding processes. 

Moreover, it was shown that under well defined processing conditions the secondary ion-beam parameters (beam 
profile) may lead to different patterns [23] • Theoretical predictions including statistical properties of the beam 
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profile were performed in Ref. [9]. It was shown that fluctuations in incident angles result in stochastic description 
of the ripple formation with multiplicative noise. Unfortunately, detailed description of pattern formation in such 
complicated stochastic systems was not discussed. Moreover, the problem of understanding the scaling behavior of 
the surface characteristics is still opened. 

In this article we aim to study ripple (or generally pattern) formation processes in anisotropic system gov- 
erned by the corresponding Kuramoto-Sivashinsky equation which takes into account multiplicative noise caused 
by fluctuation of the incidence angle. We consider the linear and nonlinear models separately and discuss the 
corresponding phase diagrams in the space of main beam parameters reduced to the penetration depth and the 
incidence angle. Moreover, we present results of the scaling behavior study of the correlation functions and discuss 
time dependencies of the roughness and growth exponents during the system evolution as well as fractal properties 
of the surface. It will be shown that multiplicative fluctuations in ripple formation processes can accelerate/delay 
surface modulations. We shall show that both phase diagrams and the scaling exponents crucially depend on the 
statistical properties of the beam. 

The work is organized as follows. In Section[2]we present the stochastic model with multiplicative noise. Section[3] 
is devoted to the stability analysis of the linear system, where the main phase diagrams are discussed. The nonlinear 
stochastic model is studied in Section |4] Here we consider the behavior of the main statistical characteristics of 
the surface such as distribution of the height field, scaling properties of the correlation functions. We summarize in 
Section [5j 

2 Model 

Let us consider a (/-dimensional substrate and denote with r the d-dimensional vector locating a point on it. The 
surface is described at each time t by the height z = h(r,t). If we assume that the surface morphology is changed 
while ion sputtering, then we can use the model for the surface growth proposed by Bradley and Harper and 
further developed by Cuerno and Barabasi [4] . We consider the system where the direction of the ion beam lies in 
x — z plane at an angle 9 from the normal of the uneroded surface. Following the standard approach one assumes 
that an averaged energy deposited at the surface (let say point O) due to the ion arriving at the point P in the solid 
follows the Gaussian distribution [3] E(r) = (e/(2n) 3 ^ 2 afj, 2 ) exp(— z 2 /2a 2 — (x 2 + y 2 )/2p, 2 ); e denotes the kinetic 
energy of the arriving ion, a and p are the widths of the distribution in directions parallel and perpendicular to the 
incoming beam. Parameters a and p depend on the target material and can vary with physical properties of the 
target and incident energy. We consider the simplest case when a — p. The erosion velocity at the surface point 
O is described by the formula v = p J„ dr$(r)_E(r), where integration is provided over the range of the energy 
distribution of all ions; here $(r) and p are corrections for the local slope dependence of the uniform flux J and 
proportionality constant, respectively |24j . The general expression for the local flux for surfaces with non-zero local 

curvature is [5Sj: <fr(x,y,h) = J cos (arctan y/JV x h) 2 + (V y h) 2 J. Hence, the dynamics of the surface height is 

defined by the relation d t h ~ — v{6 — V ffi /i, V^/i, V?ft) and is given by the equation d t h ~ —v(9)y/l + (V/i) 2 , where 
< 9 < tt/2 p3 [S7J SI 119 [SJ. The linear term expansion gives d t h = -v + jV x h + v a \I 2 aa h\ where V = d/dr, 
V Q = d/da, a = {x, y}. Here vq is the surface erosion velocity; 7 = 7(6*) is a constant that describes the slope 
depending erosion; v a — v a {9) is effective surface tension generated by erosion process in a direction. 

If one assumes that the surface current is driven by differences in chemical potential p, then the evolution 
equation for the field h should take into account the term —V • j s in the right hand side, where j s = KV(V 2 h) 
is the surface current; K > is the temperature dependent surface diffusion constant. If the surface diffusion is 
thermally activated, then we have K = D s Kp/n 2 T, where D e = Doe~ E "/ T is the surface self-diffusivity (E a is the 
activation energy for surface diffusion), k is the surface free energy, p is the areal density of diffusing atoms, n is 
the number of atoms per unit volume in the amorphous solid. This term in the dynamical equation for h is relevant 
in high temperature limit which will be studied below. 

Quantities vq, 7, v a are functions of the angle 9 only, not the temperature. Assuming that the surface varies 
smoothly, next we neglect spatial derivatives of the height h of third and higher orders in the slope expansion. Taking 
into account nonlinear terms in the slope expansion of the surface height dynamics, we arrive at the equation for 
the quantity h' = h + vot of the form (3j 0] 

d t h = jV x h + v a V 2 aa h + ^(V Q /i) 2 - ifV 2 (V 2 ft), (1) 



where we drop the prime for convenience. Coefficients in Eq.pJ) are defined in Ref.(3] and read 

s = sm9, c = cos 9, a a = a/a, F = (epJ/\f2ir) exp(— a%/2), 

7 =f s (a 2 CT c 2 -l), 

X x = £ c(al (3s 2 - c 2 ) - a\s 2 c 2 ), A, = -f c(a 2 c 2 ), 



v x = fa CT (2s 2 -c 2 -a 2 s 2 c 2 ), i/„ = -f a^c 2 



Here all control parameters are defined through the ion penetrate distance a, the incidence angle 0, the flux J and 
the kinetic energy e. It is known 25 that the penetration depth depends on the target material properties and the 
incoming ion energy e: a « e 2m /nC m , where n is the target atom density, C m is the constant depending on the 
interatomic interaction potential [2lT , to sa 1/2 for intermediate energies (from 1 to 100 keV). Equation ([Tbis known 
as the noiseless anisotropic Kuramoto-Sivashinsky equation 22 

It was shown [3] that the linearized dynamical equation (nj) admits a solution of the form h(x, y, t) = Aexp(i(k x x+ 
k y y — tot) — rt), where cj = —j(9)k x is the frequency, r = — [v x {9)k 2 + v y {9)ky) — K(k 2 + k 2 ) 2 is the parameter 
responsible for a stability of the solution. During the system evolution a selection of wave-numbers responsible for 
ripple orientation occurs. The selected wave-number is fc 2 = \i/ a \/2K, where a refers to the direction (x or y) along 
which the associated v a has smaller value. 

For the noiseless nonlinear model (TlJ) it was shown that due to the sets v a and X a are the functions of the 
incidence angle 9 £ [0,7r/2] there are three domains in the phase diagram (a a ,9) where v x and X x changes their 
signs, separately [4]. It results in ripples formation in different direction x or y varying a„ or 9. 

To describe an evolution of the surface in more realistic conditions one should take into account that the 
bombarding ions reach the surface stochastically, i.e. at random position and time; generally, it can reach the 
surface at random angle lying in the vicinity of the incidence angle 9. Most of models proposed to describe 
ripple formation due to the ion sputtering process incorporates additive fluctuations £(r, £) that takes into account 
stochastic nature of arriving ions (see for example Refs.|4, 6, 14 ). From the mathematical viewpoint such stochastic 
source results in spreading the patterns and makes possible statistical description of the system. If this term is 
assumed as a Gausisan white noise in time and space it can not change the system behavior crucially [281 129) . 

If one supposes that the ion beam is composed of ions distributed with different incidence angles, then we 
have three possible cases [3]: (i) homogeneous beam when the erosion velocity depends upon random ion beam 
parameters and the average velocity is defined through the distribution function over beam directions; (ii) temporally 
fluctuating homogenous beam when the direction of illumination constitutes a stationary, temporally homogeneous 
stochastic process; (iii) spatio-temporally fluctuating beam when the directions of ions form a homogeneous and 
stationary field. In Ref.[5] authors consider the case (iii) under assumption of the Gaussian distribution of a 
beam profile centered at a fixed angle 9$. Such model means that the fluctuation term that can appear in the 
dynamical equation for the field h is some kind of a multiplicative noise (with intensity depending on the field 
h). Unfortunately only general perspectives were reported for the nonlinear model, whilst main results relate to 
studying the linear model behavior. From the naive consideration one can expect that the multiplicative noise can 
qualitatively influence on the dynamics of ripple formation in the nonlinear system. 

In present article we aim to consider the general problem of the ripple formation under assumption of Gaussian 
distribution of the beam profile around 9q in the framework of the model given by Eq.|fij| following the approach 
proposed in Ref.[9]. To describe the model we start from Eq.Q which can be rewritten in the form dth = f(6, V a /i), 
where / is a deterministic force. Considering small deviations from the fixed angle 9q we can expand the function 
f(9, V Q /i) in the vicinity of 9q. Therefore, for the right hand side we get / = /o + (df/d9)\g=g 89 and assume that 
59 is a stochastic field, i.e. 59 = 59(r,t). Assuming Gaussian properties for the stochastic component 59, we set 

(59{r, t)) = 0, (59(r, t)59(r', t')) = 2£>SC r (r - r')C t (t - if), (2) 

where D is the parameter depending on the beam characteristics such as J, e, p, a, a; S is the noise intensity 
characterizing dispersion of 59; C r and Ct are spatial and temporal correlation functions of the noise 59. In further 
consideration we assume that 59 is the quasi- white noise in time with Ct(t — if) — > 5(t — t') and colored in space, 
i.e. C r (r — r') = (y / 27rr 2 )~ d exp(— (r — r') 2 /2r 2 ), where r c is the correlation radius of fluctuations. At S = no 
fluctuations in the beam directions (incidence angle) are realized (pure deterministic case). 

Therefore, expanding coefficients at spatial derivatives in Eq. (nl we arrive at the Langevin equation of the form 



d t h = ToV^ + ^oVL/i + ~ (V a h) 2 - KV 2 {V 2 h) 



TiV^/i + ^iVL/i + ^Vc^) 2 



(3) 



where 70 = 7(^0), ^qo = v a (0o), A q0 = A Q (0 O ), 71 = d e ^f\e=e , v<xi = d e v a \e = e 01 X al = a e A a | e= e . The parameter 
D is reduced to the constant F, that means that multiplicative fluctuations appears only if the system is subjected 
to ion beam with F^O. Therefore, the stochastic system is described by the anisotropic Kuramoto-Sivashinsky 
equation with the multiplicative noise. 



3 Stability analysis of the linear model 

It is known that transitions between two macroscopic phases in a given system occur due to the loss of stability 
of the state for the certain values of the control parameters. In the case of stochastic systems the liner stability 
analysis needs to be performed on a statistical moment of the perturbed state. We will now perform the stability 
analysis for the system with multiplicative fluctuations. To that end we average the Langevin equation (Hi) over 
noise and obtain 



d t (h) =TtoV„</») + ^oVL W + ^f ((V Q ^) 2 ) - KV 2 (V 2 (h)) 

The last term can be calculated using the Novikov theorem [30] . From a formal representation one has 
(nS9(x, y;t))=J At' J Ax' J Ay' (69(x, y; t)S0(x', y'; t')) ^ ^fry.^)) ) ' 



(4) 



(5) 



where 1Z is the functional, 5/5(59) is the variational derivative. The integration is carried out over the whole range 
of x', y' and t'. For our model one has TZ — jiVxh + UaiV^h + ^p(V Q /i) 2 . The variational derivative can be 
computed with the help of the relation j^j|y = ^ (&g) _ ,, where the second term is obtained from the formal 
solution of the Langevin equation Q. It follows that the response function takes the form 



(lite) = 7lV * M(a; " xl) + s{a - a ' ] { v «tfLh + ^f^ a hf 



(6) 



As a result the variational derivative can be written as follows 
5K 



5(59) 



=7iVx 



yi7 a h5(x - x) + 5(a - a') <j v al V 2 aa h + -f(V a hf 



Aal 



-v a \J ac 



~f\V x h8{x - x') 



r,\ m vl P ,h+^(v l3 h) 



(7) 



+X al (V a h)V a 



7iV x W(a! - x') + 5(/3 - f3') L 01 W 2 00 ,h + ^f(V phf 



Let us consider the stability of the linear system. From the relation obtained it follows that terms with coefficients 
A a i lead to the nonlinear contribution, and hence can be neglected at this stage. Therefore, reduced expression is 
of the form 



511 
5(59) 



=1 2 V X [V x h5(x 



- iiv al {V x [V 2 aa h5(a - a')] + V 2 aa [V x hS(x - x')} } 



-^i^iVL [V^M(a -/?')] ■ 



(8) 



To perform next calculations we assume that the spatial correlation function for fluctuations can be decomposed 
as C r (r) = C x (x)C y (y) with maximum at a = a', where C(0) = C x (0) = C y (0) and C"\ a=a < = d xx C x \ x=x > = 
d 2 y Cy\y =y ' , with C"\ a=a ' < 0. Then, integrating over t', and x' and y' (by parts), we obtain the expression for the 
decomposed correlator: 



< [7iV x ft + v a ^ 2 aa h +(\ al /2)(V a h) 2 ]] 59) ~ 

{^iC"| a=a 'VL + 7i 2 c(o)vL + c(o)KiVL) 2 

+71 v^ [C"| a=a /V x + C(0)VLJ +7i^iC(0)V 2 Q V ;E } (h). 



(9) 



Finally, we can rewrite the linearized evolution equation for the average (h) in the standard form: 

dt(h) = K f (h) + iZf{h)-iQ(h), (io) 

where the following notations are used 

tT/ = (7o + 71 S [u xl C"\ a= a> + ^xiC(0)VL + ^iC(0)VL] )v B , 

£) = (^0 + EC"| a=Q /i^)VL + 7i 2 SC(0)VL, (11) 

JC/ = -i^(vL) 2 + EC(o)(i/ a iVL) 2 - 



It is easy to see that Eq.(10) admits a solution of the form (h) = Aexp(i(k x x + k y y — ujt) +rt). Indeed, substituting 



it into Eq.(10) and separating real and imaginary parts we found 



w = -fc x (7o + ^\v xl Y>C"\ x=X ') + k x r yiv x i'EC(0)(k%, + k 2 y ), 

r = -k 2 x T x - k 2 y T y - K(k 2 x + k 2 ) 2 + £C(0)(i&*£ + z^fc 2 ) 2 , (12) 

r, = v x0 + v 2 xX ^c"\ x=x , + 7 2 ec(o), r y = Vya + v 2 yl Y,c"\ y=y ,. 

It follows that if T Q < 0, then there will be a range of low frequencies that will grow exponentially. From our model 
one can see that as v y o < and C" < with C(0) > the quantity T y is always negative. Therefore, instability 
along y direction will always appear. The quantity T x can change sign as control parameters 9 and a a and noise 
characteristics vary. It means that instability in x direction can appear at same incidence angles and penetration 
depths. Moreover, the statistical characteristics of the noise reduced to the spatial correlation length r c and the 
intensity S governing the total stability of the solution can change the system behavior drastically. 

Stability change of the anisotropic system with an additive noise was discussed earlier [4]. Let us consider 
stability change in the system with the multiplicative noise. In Figures [l^,b we plot the corresponding phase 
diagrams at fixed noise intensity E and different correlation radius r c . Here dotted lines limit domains of the 
stability of the system at low frequencies and relate to the case T x = 0. Solid line divides the space of a a and 9 
where parameter B = 2K — 4EC(0)f„ 1 takes zero values at k x = k y . This parameter is responsible for the stability 
of the system at large wave-numbers. It is known that observable/selected ripples correspond to wave-numbers with 
^a = |r|/-B where B > and V — mxa^xjTy]. Dashed lines in Figfll correspond to the system parameters where 
k x — k y . In domains denoted with the corresponding wave-number k x or k y ripples have the orientation in x or 
in y direction, respectively. As it follows from our linear stability analysis, orientation of ripples can be controlled 
varying both the penetration depth a a and the angle of incidence 9 at fixed E and r c . Comparing plots in Fig Ilk 
and in Fig{TJ3 one can see that the statistical properties of the noise 89 are responsible for the change of the system 
behavior. Indeed, at small correlation radius of the angle fluctuations the domain of the system instability at fixed 
a a = 1 is bigger than at large r c . Moreover, the variation of quantity r c can lead to a decrease of the system 
parameters where ripples oriented along k x are observed. It is interesting to note that at large r c at fixed interval of 
the incidence angles 9 a reorientation of ripples can be found varying parameter o CT related to the deposited energy 
of the beam. Indeed, in the interval of 9 lying between the abscissa of point E and abscissa of point F some kind of 
reentrance is observable: at small a a (below the bottom dashed line where k x — k y ) the ripples are oriented along 
k y ; in the intermediate domain of a a (between two dashed lines) the ripples are oriented along k x ; at large a a the 
ripples are oriented along k y again (see snapshots for points A — D). The same situation is realized at fixed a a when 
the incidence angle varies. For the system parameters related to the dashed lines (see points E, F) the ripples are 
characterized by k y = k y with the orientation angle 7r/4. 

Next, we calculate the selected wave-lengths A^ and X y versus the angle of incidence 9 and the penetration depth 
a a (FigJ2k) and versus the correlation scale r c and the energetic parameter F (FigJ2b). The selected wave-lengths 
relate to the smallest wave-number in the corresponding direction. It is seen that as a a or 9 varies transformations 
in ripple orientation occur. Here a a i and 9i are threshold magnitudes for the penetration depth and incidence 
angle, respectively, indicating change of the ripple orientation. It is seen that there are two critical values a a \ = a c ax 
and (Zo-5 = a° where the corresponding wave-lengths take infinitely large magnitudes due to T x — 0. There are 
two critical value for the angle 9-2 = 9 X and #3 = 9y indicating divergence of the wave-lengths when T x takes 
zero values. From Fig(2J3 one can see that as the energetic parameter F increases the wave-length of the ripple 
formation reduces to zero. At small a a orientation of selected ripples can be changed at F — F\, whereas at large 
values for the penetration depth no change is possible in the ripple orientation. The dependencies A a (r c ) manifest 
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Figure 1: Phase diagrams for the pattern selection at in the system with multiplicative noise with E = 1 (in domains 
denoted with k a where a £ {x,y\ patterns with wave- vector k = |fc Q |a are selected; plots (a), (b) correspond to 
r c = 0.65, r c — 1, respectively). 



non-monotonic behavior: at small r c the wave-length increases, whereas at large r c the decreasing dependencies are 
observed. Moreover, there is a critical value for the correlation radius r cl where orientation of ripples can be changed. 
Therefore, correlation properties of the ion beam can play a crucial role in ripple formation processes at early stages 
(in linear models). From the equations obtained for the selected wave-numbers it follows that the selected wave- 



lengths have the well-known assymptotics versus main parameters of the beam (A ' 
A ~ J" 1 ' 2 ) and depend assymptoticaly versus secondary characteristics: A ~ (Sq 



T-Va exp( 



-E a /T),\~e 
(?'c0 -Tc^ 1 . 



-1/2 



4 Nonlinear stochastic model 

Next, let us consider the nonlinear system behavior setting A Q ^ 0. In further study we are based on the simulation 
procedure allowing us to solve the nonlinear stochastic differential equation Q. As it was done in previous section 
we use the finite-difference approach to calculate the evolution of the field h. 



4. 1 Evolution of the height distribution function 

To investigate properties of a distribution of the field h we use skewness m^ and kurtosis mi, defined as 

„ {(h(r)-{h(r))) 3 ) 



»*3 — 


w 3 

((Mr) - (Mr))) 4 ) 


w 2 -- 


= ((h(r)-(h(r)}) 2 ), 



(13) 



where (h(r)) is the average of the height field ((/i(r)) = V^ 1 J2 r M r ^)> V = L d is the system volume, d is the 
spatial dimension, L is the linear size of the system); W is the interface width. Skewness is a measure of the 
symmetry of a profile about the reference surface level. Its sign tells whether the father points are proportionately 
above (7713 > 0) or below (7713 < 0) the average surface level. Kurtosis describes randomness of the surface relative 
to that of a perfectly random (Gaussian) surface, for the Gaussian distribution one has 777,4 = 3.0. Kurtosis is a 
measure of the sharpness of the height distribution function. It is known that if most of the surface features are 
concentrated near the mean surface level, then the kurtosis will be less than if the height distribution contained a 
larger portion of the surface features lying farther from the mean surface level. 
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Figure 2: (a) Plot of selected wave-lengths X x and Aj, vs. incidence angle at a a 
penetration depth a a at 9 = 0.4763, r c — 1.0 (a a i with i e 1, . . . ,4 denotes threshold values when a change of the 
wave- vector of patterns occurs); here < 9 < 7r/2 is measured in radians; other parameters are: F = 1, cr = 1. (b) 
Plot of selected wave-lengths X x and A,, vs. the correlation scale r c at F = 1.0 and the energetic parameter F at 
r c = 1.0 at 9 — 0.4763. Other parameters are: F = 1, a = 1, £ = 1. 



Figure [3^, shows snapshots of the surface morphology for the set of parameters: a a — 1.2, 9 = 0.4, F — 1.0, 
a = 1.0, K = 2.0, S = 1.0, r c = 1.0 at t = 20, 40, 60, and 100, respectively. In our simulations we have used 
Gaussian initial conditions taking (h(r,t = 0)) = 0, ((5h) 2 ) — 0.1; integration time step is At = 0.005, space step is 
I = 1. It is seen that with an increase of the growth time, the lateral length of the surface features becomes bigger 
and holes (black regions) are formed at t = 100. It follows that due to nonlinear effects and noise action the surface 
morphology is crucially changed comparing to the ripples shown in Figllb. Figure [3b illustrates the corresponding 
height probability density distribution function of these surfaces. It is seen that at t = 20 the distribution is close 
to the Gaussian distribution. With the increase of the growth time, there is deviation from zero-centered Gaussian 
distribution and after transient period of time the probability density function becomes symmetrical and centered 
around zero. In Figpt we plot the kurtosis m^, the skewness m 3 , and the interface width W as functions of the 
growth time for above system parameters. According to initial conditions we have m^ ~ 3.0, m^ ~ and W ~ at 
t ~ 0. With the increase of the growth time the kurtosis grows until maximum is reached. The skewness decreases 
to its minimum, and after tends to zero. These two quantities reflect the form of the distribution function shown 
in Fig(3f). The interface width increases algebraically toward a saturation regime at large t. 

We have computed phase diagram for the nonlinear systems illustrating formation of different patterns shown 
in Fig|4j It is seen that the numerical results are well related to analytical predictions from the linear stability 
analysis. Indeed, critical points lying on the lines correspond to a change of the sing of the quantity r^. At large 
and low penetration depth a a ripples oriented along k y direction are observed (see snapshot for the point C) due 
to r^ > 0. At the intermediate values of a a random patterns (holes) are realized due to the nonlinear influence of 
both the deterministic term X a (V a h) 2 and the stochastic contribution. 



4.2 Scaling properties of the surface morphology 

Using numerical data it is possible to study statistical properties of the system considering the time-dependent 
height-height correlation function, determined as follows C/j(r,t) = ((h(r + r',t) — h(r',t)) 2 ). In the framework of 
dynamic scaling hypothesis one can write the correlation function in the form |31l I32j 



C h (v,t) = 2W 2 (t)<i> 



m 



where 




for ii«l, 
const, for m>1. 



(14) 
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Figure 3: A typical evolution of the system with multiplicative fluctuations at early stages, (a) Snapshots of images 
of the field h distribution for various growth times (dark color indicate low h, white areas relate to high h). b) 
Probability density function of the height for various growth time, c) Kurtosis, skewness, and interface width versus 
growth time. Other parameters are: a„ = 1.2, = 0.4, F = 1.0, a = 1.0, K = 2.0, £ = 1.0, r c = 1.0. 



Early stages can be fitted by the function [33J CVi(r,£) ss 2W 2 (t)[l — exp[— (r/£,) 2a ]. The dynamic scaling hypothesis 
assumes that the following dependencies are hold: W 2 (t) ex t 2/3 , £(£) oc t 1 / 2 , where /3 is the growth exponent, z is 
the dynamic exponent for which z = a/ j3. From another viewpoint one can assume |34j 



C h (r,t)=r 2a i^ 



(16) 



where 



■ ,2rf 



for u«l, 



m~< ' . ::: (17) 

I const, for v >• 1, 

and the relation z = a//3 holds. Therefore, these two cases lead to the same results Ch(t) oc t 2 ^ and Ch{r) oc r 2a , 
allowing one to define the growth exponent /3 and the roughness exponent a. As was shown in Ref . [34) the 
roughness W(t,L) can be related to the structure function S(k) as follows W 2 (t,L) = V^ 1 ^ k , Q S(k, £), where 
Sh(k,t) = V^ 1 (h k (t)h_ k (t)) . The structure function S(k,t) has the form 

S h (k,t) = k- {d+2a ^G(k z t), 



where 



@(k z t) 



k 2a t 2a /P, for fc z t< 1, 
const, for k z t 3> 1, 



(18) 



(19) 
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Figure 4: Phase diagram for the anisotropic nonlinear model at F — 1, a — 1, 'S = 1, r c = 1. Snapshots are taken 
at the system parameters related to position of the points A, B and C, respectively. 



and scales as Sh(k,t) ex fc~( d + 2Q ) for large t and Sh(k,t) ex t 213 for small t. 

In previous studies (see, for example Ref. :6) it was shown that even in the isotropic system with additive noise 
scaling exponents a, j3 and z depend on the system parameters vq, Ao and K. Moreover, these exponents are the 
time-dependent functions, i.e. its magnitudes can be changed in the course of the system evolution. 

In our study we have taken into account multiplicative noise described by the energetic characteristics of the 
beam and additionally by the noise intensity £ and correlation radius of fluctuations r c . Therefore, one should 
await that due to renormalization of the main system parameters responsible for the stability of the system and 
nonlinear effects in its behavior such scaling exponents are functions of the above noise properties. To prove it we 
compare magnitudes of both scaling exponents a and j3 for the system with additive fluctuations and for the system 
with our multiplicative noise. 

According to the scaling hypothesis the temporal evolution of the quantity W = ((5h) 2 ), where 8h = h— (h), can 
be represented through the exponent 8 related to the exponent f$ as 6 = j3. It is known that the ordinary diffusion 
(Brownian) process is described by Einshtein law ((Sh) 2 ) ex t 2S , with 8 = 1/2. If the exponent 8 deviates from the 
value 1/2, then anomalous processes are realized: at < 8 < 1/2 there is a delayed (subdiffusion) process, whereas 
at 8 > 1/2 the accelerated diffusion (superdiffusion) is realized. By comparison of results related to additive and 
multiplicative noise influence in anisotropic system, one can see that in the case of the additive noise influence 
we get 1/2 < 8 < 1. In the case of the multiplicative noise influence the quantity 8 takes values in the window 
< 8 < 4. It means that at small time intervals there are delayed processes which can be accelerated by the noise 
action at intermediate t; at large t a transition toward saturation regime reduces growth velocity, decreasing 8. 

To characterize fractal properties of the surface one can study a pair correlation function defined as follows: 



C v (v;t) = (h{r + r';t)h(r;t)). 



(20) 



If there is no characteristic space scale, then the introduced correlation function should behave itself algebraically, 
i.e., C p (r;t) ex l/r A , where the scaling exponent A relates to the fractal correlation dimension Di as A = d — Z?2- 
The corresponding Fourier transformation of the correlation function C p (r;t) scales as S p (k;t) ex k~ D2 . From the 
definition of the correlation fractal dimension £>2 and the properties of the Fourier component of the correlator 
(20) it follows that at D-i = there is no scaling behavior of the structure function and S p (k;t) sa const. Hence, 
the surface at the fixed time t can be considered as a Gaussian surface with no correlation, i.e. white noise in 



space with equal contribution of all wave-numbers k; the corresponding spatial correlator (20) is reduced to the 



Dirac delta-function, C p (r) — > S(r). In the case D-i = 2 one arrives at typical dependence S p (k) ex k 2 for diffuse 
spreading on the structured (let us say, flat) surface. Here the topological dimension d equals the fractal dimension 
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Figure 5: Roughness exponent a and growth exponent /? versus growth time for isotropic Kuramoto-Sivashinsky 
equation with additive noise (white circles and squares) at v x = v y = —0.2, X x = X y — 1.0, K — 2 and anisotropic 
Kuramoto-Sivashinsky equation with additive noise (black circles and squares) at a a = 1.2, 9 = 0.4 K = 2. 
Snapshots are shown for above two models (left for isotropic and right for anisotropic Kuramoto-Sivashinsky equation 
with additive noise) at t = 20, 60, 100 from top to bottom. The noise intensity £ = 1.0. 



D 2 • Therefore, a variation of the fractal dimension D 2 versus the time indicates a change of the fractal morphology 
of the surface from pure uncorrelated Gaussian surface toward well structured surface having fractal dimension 
d = D 2 = 2. 

In order to study scaling properties of the system under consideration we will compare all our results obtained 
with results coming from an investigation of the anisotropic system with additive fluctuations. Such system will 
serve as a reference system. Moreover, to verify and to test our the numeric procedure of the scaling exponents 
computation we recalculate results of the work [5] for the isotropic Kuramoto-Sivashinsky equation. 

As a reference system we consider the model described by the Langevin equation with additive noise, i.e., 
dth = i/ a oV^ ta h+ ^f-(V a h) 2 — KX7 2 (X7 2 h) + C(r, t), where £ is the Gaussian random source with properties (£) = 0, 
(C(r, t)((r' , t 1 )) = 2£<S(r — r')S(t — t'). Calculations of the dynamical exponents at the system parameters a a = 1.2, 
9 = 0.4, F = 1.0, a = 1.0, K = 2.0, £ = 1.0, r c = 1.0 are shown in Fig(5] It is seen that the exponents a and (3 for 
above two models are different. In the anisotropic case we have elevated magnitudes for a and /3, i.e. such exponents 
essentially depend on the control parameters of the system. Hence, due to renormalization of the control parameters 
by the multiplicative noise contribution the dynamic scaling exponents depend on the noise characteristics. 

Let us consider the anisotropic system with multiplicative fluctuations. We have performed calculations of the 
scaling exponents at the fixed point on the phase diagram (9,a a ) at different values of the noise intensity S and the 
correlation radius r c . The reference point is a a = 1.2, 9 = 0.4, F = 1.0, a — 1.0, K — 2.0. We compute a and /? at 
time window when the interface width W or the correlation function Ch(r) start to grow until they saturate (i.e., 
when algebraic dependencies W 2 (t) ex t 2/3 and Ch(r) ex r 2a are observed). 

The corresponding time dependencies of a, (3 and D 2 are shown in Fig(6j In Fig(6^, we plot the corresponding 
correlation function C/j(r; t) and the roughness exponent a; Figurepp illustrates the time dependence of the interface 
width W and the growth exponent (3; Figure|6t shows the pair correlation function C p (r; t) and the associated fractal 
dimension D 2 at fixed times. From Figj6^, it is seen that the growth process is nonstationary for early stages and 
the roughness exponent a is near 0.95 for small incidence angle dispersion E. At such set of the control parameters 
(a a and 9) the correlation radius r c has not essential influences on the system behavior. At large £ the roughness 
exponent has lower magnitudes and a has the well pronounced time dependence. 

Comparing curves for the interface width at different £ and r c from the one hand and the growth exponents 
dependencies versus time from another one (see FigJ6p) , one can conclude that as the noise intensity £ increases the 
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Figure 6: (a) The correlation function Ch(r) and the roughness exponent a versus time, (b) The interface width 
W and the growth exponent f3 versus time, (c) The correlation function C p (r) at t = 60 and the fractal dimension 
D2 versus time. Other parameters are: a a — 1.2, 9 — 0.4, F = 1, a = 1, K = 2. 
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position of the peak of the exponent /3 reduces to small time. It means that as the noise intensity increases at such 
choice of the control parameters the interface width increases at smaller time interval than at low E. Alternatively, 
the shift of the peak position at large E indicates that multiplicative fluctuations are responsible for nonlinear 
effects at small times. It looks natural due to the nonlinear form of the multiplicative noise, where the large noise 
contribution influences crucially on properties of the growth processes. The correlation properties of fluctuations 
characterized by r c define a height of the peak for /3. In other words, the noise correlations can accelerate growth 
processes increasing the interface width W until it attains the saturation. 

A change of the fractal properties of the surface is shown in Fig(6b. Here we compare fractal properties of two 
different systems, namely with additive fluctuation source and with multiplicative noise. From the dependencies of 
the pair correlation function C p (r) it is seen that at t = 60 the additive noise contribution leads to a picture when the 
correlation function C p (r) decreases slowly with exponent A = 0.227, whereas the multiplicative noise contribution 
with the same intensity E = 1.0 at r c = 1.0 increases the exponent A to 1.587. According to the definition of the 
correlation dimension D 2 it means that the fractal properties of the surface is well pronounced at multiplicative 
noise with large intensity at a small time interval t ~ 60 (see curves D%{t)). From the time dependencies of the 
fractal dimension D 2 for the system with multiplicative noise it follows that at small times the surface has Gaussian 
properties of the kind of white noise in space (the correlation function has the from of the Dirac delta-function.). 
At small time interval (at intermediate times) the fractal properties emerge and characterized by < D 2 < 2. At 
large times one has D<i = 2 and the well structured patterns are observed, its dimension D 2 coincides with the 
topological d — 2. In the case of additive fluctuations the time interval of the formation of well structured patterns 
is larger than in system with multiplicative noise. 

Next let us compare the time dependencies for the scaling exponents for different set of the system parameters 
a a and 9 shown in Figs[7k,b. It is seen that at a a — 2.0, 9 = 0.4 (FigjTk) at large noise intensity E the interface 
width grows faster comparing with the case of small E at fixed r c . But the growth at small E occurs at earlier 
periods of time. This situation absolutely different to the case shown in Fig|6b, whereas variations in r c leads to 
the increase of the maximum for j3. The roughness exponent a does not principally change its values at different 
t. Comparing curves for the fractal dimension D 2 for above two sets of the system parameters one can see that 
at small E the quantity D 2 is smaller than in previous case (cf. dependencies D 2 (t) in Figspfc, [7k). The same 
dependencies of the scaling exponent a, (3 and D 2 at a„ = 2.0, 9 = 0.2 are observed. Here at large E the nonlinear 
effects delayed. Considering snapshots shown in the right hand side of the plots at different E and r c one can see 
that depending on the noise properties the morphology of the surface changes crucially (cf. snapshots at different 
r c in FigJ7p). It means that the phase diagram shown in FigHwill be modified under variation of E or r c as linear 
stability analysis predicts. 

Figure [8] show the evolution of the spherically averaged structure function defined on a circle 

S h (k,t) = ±- Yl ^(k,i), (21) 

fc k<k<k+Ak 

where iVj. is the number of point on the circle of the width Afc. Two different choice of the noise intensity E and 
the correlation radius of fluctuations r c are shown in Figs(8^,,b,c, respectively. It is seen that during the system 
evolution at early stages the system selects the ripples with the corresponding wave-number (dotted lines) and after 
at late stages the algebraic form for the structure function is realized. At large time intervals one can define the 
exponent a from the definition Sh(k) ex k~( d+2a \ As is shown in Figsl8a,b for above system parameters one has 
Sh(k) oc fc~ 3 - 86 for E = 0.1 and Sh(k) ex fc~ 3,92 for E = 1.0, where d = 2. Hence, the roughness exponent takes 
values a ~ 0.93 and a ~ 0.96 that is well predicted by the analysis of the correlation function Ch(r) (see FigsJ6fe 



5 Conclusions 

We have studied the ripple formation processes induced by the ion sputtering under stochastic conditions of illu- 
mination. The main assumption was stochastic nature of the ion beam when the angle of incidence distributed 
in space and time (homogeneous and stationary field). It allows us to generalize the Bradley-Harper model of 
ripple formation [3] and consider the stochastic model with multiplicative fluctuations describing random nature of 
the incidence angle proposed in Ref.^S]. We have discussed properties of the ripple formation in both linear and 
nonlinear models. 
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Figure 7: Scaling exponents a, /3 and correlation dimension D2 versus time at (a) a a = 2.0, 9 — 0.4 and (b) 
a a = 2.0, = 0.2. Typical snapshots of the system evolution are taken at different S and r c . Other parameters 
are: F = 1, a = 1, K = 2 

Within the framework of the linear stability analysis we have shown that even in the linear system the noise 
action is able to change the critical values for the control parameters of the system such as the penetration depth 
and the averaged incidence angle. It was found that as correlation properties of such multiplicative noise as the 
dispersion in the incidence angles around the average can reduce the domains of the control parameters where the 
ripples change their orientation at the fixed angle of incidence. 

Studying the nonlinear model we have computed the dynamic phase diagram illustrating formation of different 
patterns (ripples and holes) which relates to the results from the linear stability analysis. Main properties of the 
ripple formation were studied with the help of the distribution function of the height and its averages reduced to 
the skewness, kurtosis and interface width (dispersion). To make a detailed analysis of the ripple formation we 
have examined scaling behavior of main statistical characteristics of the system reduced to the correlation functions 
and its Fourier transforms (structure functions). It was shown that as the growth and roughness exponents depend 
on the control parameters and are time-dependent functions (it was predicted by previous study of the isotropic 
Kuromoto-Sivashinsky equation [6]); these exponents depend on the noise properties: its intensity and the spatial 
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Figure 8: Evolution of the spherically averaged structure function Sh(k,t) at different noise intensity: a) £ = 0.1, 
r c = 1; b) £ = 1.0, r c = 1; c)S = 0.1, r c = 0.65. Other parameters are: a CT = 1.2, = 0.4, F = 1, a = 1. 

correlation radius. Comparing results for the system with additive and multiplicative fluctuations it was shown 
that multiplicative noise can crucially accelerate processes of ripple formation, increasing the growth exponent. 
As far as our system is anisotropic the noise action is different at different set of the main control parameter 
values. Studying fractal properties of the surface we have calculated the fractal (correlation) dimension as the 
time-dependent function. It was shown that in the system with multiplicative noise the fractal properties appear 
at small time interval of the surface growth, whereas in the system with the additive noise this time interval is 
wider. These results are well predicted by the correlation functions analysis and by Fourier transformation of the 
numerically calculated surface. 

Therefore, as patterning as the scaling behavior of the system can be controlled by additional set of parameters 
reduced to the dispersion of the angle of incidence and the correlation properties of its fluctuations. 
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